Double-heterostructure cavities: from theory to design 
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We derive a frequency-domain-based approach for radiation (FAR) from double-heterostructure 
cavity (DHC) modes. We use this to compute the quality factors and radiation patterns of DHC 
modes. The semi-analytic nature of our method enables us to provide a general relationship between 
the radiation pattern of the cavity and its geometry. We use this to provide general designs for 
ultrahigh quality factor DHCs with radiation patterns that are engineered to emit vertically. 



PACS numbers: 42.70.Qs, 42.25.Bs 
I. INTRODUCTION 

In the 25 years since their conception [HE], photonic 
crystals (PCs) have become an indispensible tool in mod- 
ern photonics experiments. One of the principal uses of 
PCs is in the creation of ultra-high quality factor (Q) 
micro-cavities, which allow the trapping of light for long 
periods of time in very small modal volumes (V) [3HS]. 
The large ratios of Q/V that can be achieved with PC 
cavity modes have enabled a range of applications requir- 
ing strong light-matter interactions, such as cavity quan- 
tum electrodynamics (QED) experiments [TUTU], optical 
switching [11], sensing [12] [13] and multiple-harmonic 
generation [T4| ITS] . 

The current state-of-the-art in PC cavities is the 
Double-Heterostructure Cavity (DHC) [4 . DHCs are im- 
portant because of their ultra-high Q factors and also 
because the cavities can be readily integrated with other 
photonic components. DHCs are realized by increasing 
the average refractive index of a Photonic Crystal Wave- 
guide (PCW) slab in a strip-like region (see Fig. []Ja)). 
This increase often involves a small perturbation, and 
is achieved by manipulating the geometry of the PCW 
lattice [4] [5] [16] , or by increasing the refractive index of 
either the slab [13 HI] or the holes [HMTl . Within the 
perturbed region the modes in the PCW are above cutoff 
and so can propagate, however in the unperturbed PCW 
the modes are evanescent, thus leading to trapping of the 
guided mode. Recent experiments have demonstrated 
DHC Q-factors of Q > 10 6 [22]. 

The development of planar PC cavities was accompa- 
nied by a discussion on the exact mechanisms of loss in 
these structures, and, by implication, how best to manip- 
ulate the Q. All planar PC cavity modes are intrinsically 
lossy because they radiate in the out-of-plane direction 
[23] ; in 2003 Noda et al. introduced the idea of gentle 
confinement, in which it was argued that the radiation 
from the cavity, and hence the loss, could be computed 
using the overlap of the Fourier components of the cavity 
mode with the light-cone. To maximize the Q, the cavity 



can be constructed to have a mode with a Gaussian-like 
envelope with a width of several periods [3] [4] [24] , in or- 
der to minimize the extent of the cavity mode in Fourier 
space. In the context of planar PC cavities this is anal- 
ogous to the work of Englund et al. [25] , and Vuckovic 
et al. [26], who showed that the far-field properties of a 
cavity mode can be computed using the fields above a PC 
slab. However, Sauvan et al. questioned the validity of 
this approach, arguing that the dominant radiative loss 
occurred due to waveguide impedance mismatch at the 
cavity boundaries, leading to Fabry- Perot reflections that 
dominate the loss [27]. The understanding of mode con- 
finement and radiation in planar PC cavities has thus-far 
been hampered by the lack of a comprehensive theory of 
confined states in these structures. 

Here, we present a first-principles theory of DHC 
modes in 3D planar PC structures, an approach that we 
designate the Frequency Approach to Radiation (FAR) 
method. In a recent short publication, we used the FAR 
to compute the Q factors and radiation patterns of DHC 
modes [28] . In this paper we provide a detailed descrip- 
tion of the FAR and use it to provide simple designs 
for DHCs with radiation patterns which have been engi- 
neered to emit predominantly in the vertical direction. 

The FAR consists of two parts: (i) we use a truncated 
basis of bound PCW modes that lie outside the light 
cone to construct a non-radiating approximation for the 
DHC mode. We use a Hamiltonian method for our mode 
expansion and generalize our previous theory [28] by in- 
cluding non-rotating wave terms, (ii) We apply pertur- 
bation theory to this non-radiating approximate DHC 
mode to compute the Fourier components of the polar- 
ization field P that lie inside the light cone. These are 
then used to compute the far-field radiation pattern and 
Q factor [29 . This method is efficient as it splits the 
task of solving a computationally intensive problem into 
solving smaller problems that are straightforward to com- 
pute. Bound PCW modes can be computed rapidly using 
well-established numerical methods [30] , and once these 
are known, the perturbation theory takes approximately 
15 minutes using our MATLAB code. This method's 
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FIG. 1: (Color online) (a) Schematic of a DHC. In the shaded 
region the average refractive index is higher than in the re- 
mainder of the PCW. (b) Schematic of a finite segment of a 
PCW that is periodic with period d in the x-direction (the 
finite out-of-plane thickness is not shown). The unit cell is 
the region between the broken lines and has period d. 



computational efficiency derives from avoiding the need 
to compute radiative modes directly. The FAR is semi- 
analytic and provides insight into the nature of the ra- 
diative processes of planar cavity structures. A central 
feature of the FAR is an integral equation that contains 
a driving term relating the geometry of the DHC to its 
modes' radiation fields, which contain much of the qual- 
itative features of the radiation pattern. Through an 
examination of this term, we provide several designs for 
realizing ultrahigh Q factor DHCs with modes that radi- 
ate predominantly in the vertical direction. 

This paper is structured as follows: Section|ll|describes 
how we expand the cavity modes in a basis of bound 
PCW modes. The formulation for computing the radi- 
ation properties of the DHC is then outlined in Section 
jHj with results presented in Sectio n |IV| and discussed 
further in Section [Vj Then in Section I VII we use the FAR 
to design DHCs which radiate predominantly in the ver- 
tical. In Section [VlII we discuss our results and conclude. 
Appendix [A] outlines our approach for numerically solv- 
ing the integral equation which is central to FAR. 



II. NON-RADIATIVE APPROXIMATION 

In this section we describe how we obtain a non- 
radiating approximation for cavity modes. We outline 
how we expand a DHC mode using PCW modes, and 
then show that this provides a good approximation for 
the mode profile of the DHC mode. 



A. Hamiltonian formulation 

We expand the cavity mode in PCW modes using the 
Hamiltonian formalism of Sipe and co-authors [3TH33] . 
This method has proven to be useful in deriving quan- 
tum optical versions of linear [32] and nonlinear [31 cou- 
pled mode equations in a systematic way, as well as in 
devising a quantum optical treatment of dispersion and 
absorption [34]. Here, we are not interested in quantizing 
the field, but use the formulation as a vehicle for carrying 



out the field expansion and determining a first approxi- 
mation for the cavity mode. A further advantage in our 
application is that it uses the divergence-free D and B 
fields as the primary fields for the basis functions, so that 
any superposition is also divergenceless. 

We start with the macroscopic Maxwell equations 



<9D 

~dt 



V x H, 
V-D = 0, 



<9B 

~dt ~ 
V B 



-V x E, 
0, 



(1) 



with constitutive relations for non-magnetic dielectric 
media 



D(r,t) =eoE(r,t) + P(r,t), 
B(r,t) =MoH(r,i). 



(2) 



We first consider the unperturbed structure, i.e. a dis- 
persionless, non-absorbing PCW with an isotropic linear 
response. Since the DHC geometry involves perturbing 
the underlying PCW (Fig. []Ja)), the PCW modes form 
a natural basis to expand the DHC modes. Taking D to 
be one of the primary fields, the polarization field is 



P(r,t)=r(r)D(r,t), 



(3) 



where f (r) = (e(r) — l)/e(r) and e(r) is the permittivity 
distribution of the PCW. The Hamiltonian representing 
the energy of the field is 



H 



1 

2/lQ 



drB(r) -B(r) 



1 

2^ 



This is a canonical formulation of electromagnetism and 
commutators are introduced to produce the appropriate 
dynamics [32 . The equal time commutation relations are 



[D i (r),D^r')] = [B i (r),B^r')]=0, 
[D i (r),B^r')}=ihe^ k ^[d(r-r') 



(5) 



where the superscripts ij,k indicate cartesian compo- 
nents, e^ h is the permutation symbol, and repeated su- 
perscripts indicate summation. The dynamics of the 
fields are given by the Heisenberg equations of motion 



ih 
ih 



dt 
dB 

dt 



[B,H]. 



(6) 



Equations Q-Q reproduce the Maxwell curl equations 
in 0. The divergence equations in ([!]) act as initial con- 
ditions, and if satisfied at some time, the dynamic equa- 
tions ensure that they are satisfied at all times. Note that 
in a classical framework the commutators are replaced by 
Poisson brackets (with appropriate factors of ih) and the 
operators become amplitudes. 

We take the PCW to point in the x-direction (see 
Fig. [TJb)), so the permittivity satisfies 



e(r + dx) = e(r), 



(7) 
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where d is the period of the PCW. Solutions to Maxwell's 
equations are then Bloch modes with band index m and 
Bloch wavevector k. The Bloch modes of the PCW form 
a complete set and can be used to expand any field 



D(r, t) = V / dkJ ^ia m>fc D m , fe (r)e-^™. fct + c.c, 
m J bz V 2 

(8) 

where the integration is over the Brillouin zone (BZ). 
This expansion includes all modes: those bound to the 
PCW, those bound to the slab but not the PCW, as well 
as the continuum of radiative modes that are not bound 
to the slab. The Bloch modes take the form 



D m>fe (r) 



j-d myk (r)e ikx , 



(9) 



where d m ^(r) is periodic with period d. The magnetic 
field B can be written using similar expressions to (JsJ) 
and ([9|. The normalization condition for the D field is 

dr '■ — = d mm /d(fc -k), (10) 

e e(rj 

where the integration is over all space. Using the Poisson 
summation formula, it can be seen that the field is also 
normalized over the unit cell 



, d^ >fe /(r).d m , fc (r) 
dr — = o„ 



(11) 



/ceil e oe(r) 
The normalization conditions for the B field are the same 



as in (10) and (11), but with eoe(r) replaced by /iQ. Using 
Eqs. ( 8p-( 11 ) with the commutators in Eq. ([5|, it can be 
shown after some manipulation that the operators a m ^ 
and k in (jij), where f denotes the Hermit ian conjugate, 
satisfy the commutation relations 



[ fl m,fc, flm',fe'] = 0, W m & a ln> ,k>\ = 

The Hamiltonian can then alternatively be written as 



(12) 



(13) 



which is that of a set of harmonic oscillators. 

We include the nRW terms using superpositions of the 
Bloch modes forming standing waves as basis function, 
rather than the individual Bloch modes themselves. To 
avoid double counting we discretize the BZ into an even 
number of equally spaced points N with an equal number 
of points on the positive and negative halves. This way, 
the BZ centre and edge are avoided and the points closest 
to them are k = ±ir/{Nd) and k = ±(n / d ± ir / (N d)) 
respectively. By writing the wave equation for a PCW 



V x 



V x B m?/e (r) 



e(r) 

with D m , fe (r) 



—) B m)fe (r), 

i 

-V x B m>fe (r), 



(14) 



we note that both the forward propagating Bloch mode 
B m> fc(r) and the complex conjugate of the backward 
propagating Bloch mode B^ _ fe (r) satisfy Eq. ( 14) (since 
uj m ,-k — ^m,k)- A complex conjugation of tne Bloch 
mode definition ([9| is equivalent to taking a backward 
propagating mode. This means that, apart from an over- 
all phase factor, the two modes B m> fc(r) and B^ _ fe (r) 
are equivalent. From the second of Eq. ( [l4| ), it is clear 
that this is also true for D m ^(r) and we can write 



D* 



_ fe (r) = e^D m ,,(r). 



(15) 



Re-expressing the terms in Eq. (|8| with negative k in 
terms of complex conjugated Bloch modes with positive 
k and using uj_ k = we obtain 



D(r) = W dk 

m Jk>0 



fajJr, 



-(a n 



c ) D m,/c(r) 



E 



f dk J^I {a l +a ra ,_*e- 1 *».')D^(p). 

(16) 



The superposition of two equivalent, but counter- 
propagating waves gives a standing wave which can be 
made a real function, manifested here by the addition 
of a function and its complex conjugate. We define the 
standing wave basis 



(17) 



where these functions are two orthogonal sine and cosine- 
like functions. In terms of these modes Eq. (fl6|) becomes 



C m ,fc(r) 


1 


" 1 1 




D m ,k(r) 


_ S m> fc(r) _ 


2 


—i i 




_D^(r)_ 



D(r) = E / dkV2Q c 

m Jk >° 
m Jk >° 

where we define new operators through 



(18) 



Qc,m,k 
Q s,m,k 



m,k 



1 1 

i i 

—i i 

1 -1 



Q"m,k 

-k C 
t 

m 

-ke 



a 



k 

(19)" 

As the notation suggests, Q and P act like canonical po- 
sition and momentum operators respectively. The con- 
sequences of using modes from only half of the Brillouin 
zone (k > 0) is having to use two sets of independent 
operators, i.e. those with c subscripts and those with 
s subscripts. It can be shown from Eq. (12) that these 
satisfy the commutation relations 



[Qc,m,k i Q s 



'] = 0, [P c ,, 



>] = o, 



[Qc,m,k i Ps,m' ,k'\ = 0, [Qs,m,fe)^c,m',fe'] = 0, 
[Qc,m,k: Pc,m' ,k'\ ~ ^^mm'^kk'i 
[Qs,m,k i Ps,m' ,k' ] ^^^mm'^kk'' 



(20) 
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From the new operators, potential and kinetic energy-like 
terms can be defined such that the Hamiltonian in (13) 
can be re-written as 



rn , k 



), (21) 



where p G {c, s}. This is equivalent to ( 13 ), but now writ- 
ten in terms of the new operators P and Q. In compact 
form, the mode expansion (fl8|) becomes 



(22) 



D(r,t) = 



where the index a = (p, m, fc), i.e. for brevity we use a 
sum over a to replace the sum and integral in Eq. ( [l8| ), 
and F a (r) = C m , fe (r) if p = c and F a (r) = S m ,/ C (r) if 
p = s. 

Thus far we have expressed the modes of the unper- 
turbed PCW using a new notation and we have shown 
that it is consistent with the Hamiltonian in (13). We 



now show that this formulation enables us to construct a 
Hamiltonian for the DHC while keeping the nRW terms. 
After introducing the cavity into the PCW the Hamilto- 
nian is written 7-L cav — T~L + V, where 



V=— [ dr ( - . 1 7 - - -^-) D(r) • D 
2e J \e(r) + e(r) e(r) J KJ 

= J dr 7 (r)D(r) -D(r). 



(23) 



Substituting the mode expansion (22) into Eq. (23), the 
perturbation term takes the form 



V = 2^u a ujj3Q a Qj3 J dr 7 (r)F a (r)-F c 

a, (3 



(r), (24) 



In terms of the canonical position and momentum the 
Hamiltonian is 



%cav — -Z ^2 ^ + ~ QaLafiQp, 



a, (3 



where 



L a (3 = U^SaP 



■4Lj a up J dr 7 (r)F a (r) • F^(r). 



(25) 



(26) 



The coupling of the basis PCW modes is manifested by 
the off-diagonal terms in L a p. We note that the product 
QaL a /3Q(3 includes the nRW terms. 

Double-heterostructure cavities are formed by slightly 
perturbing the refractive index of a PCW. Thus DHC 
modes have frequencies near the edge of the PCW band 
[4], and the modal fields extend over many period in 
the direction parallels to the PCW [24], so their Fourier 
transform is strongly localized around the BZ-edge. This 
means that the introduction of the cavity only couples 
PCW Bloch modes with k values near the BZ-edge. 
The frequency of DHCs being near the PCW band-edge 



also means that PCW modes of different bands couple 
weakly. Thus, to good approximation, we can expand 
DHC modes using only bound modes from the even PCW 
band. Such a band is shown in Fig.[2|a) and the field pro- 
files of some modes in this band are shown in Figs. [2jb)- 
(e). These Bloch modes are the basis functions from 
which we construct the DHC mode, and their field pat- 
terns thus control the shape of the cavity mode. As dis- 
cussed in Sect. |VII[ the variations of these Bloch are im- 
portant for the cavity mode's far-field properties. 

We compute the discrete elements of L a p in (26) nu- 
merically. If the BZ is discretized such that ofThe N 
points in the BZ, there are M/2 positive values of the 
Bloch wavevector k > below the light line, the sums 
over a contain M terms. L a p in (26) is then a real, sym- 



metric M x M matrix. Its eigenvalues and eigenvectors 
correspond, respectively, to the square of the frequency 
of the cavity modes and the associated real- valued ampli- 
tudes of the F a (r). Diagonalizing L results in a discrete 
spectrum of bound cavity modes, as well as a continuum 
of waveguide states that are not bound to the cavity. 
We are interested only in the fundamental cavity mode, 
which corresponds to the lowest eigenvalue. 

To diagonalize the eigenvalue equation, we write it as 



E 

13 



(27) 



and define an orthogonal matrix of eigenvectors 







s? ■ 


JM)1 
■ ■ *1 






4 2) . 


• • *2 


s = 












s (2) 








b M ■ 





(28) 



where S T S = SS T = I. We can therefore write 
LS = SO, 



(29) 



where Vt = diag(cl^) is a diagonal matrix of cavity mode 



frequencies. Defining new canonical coordinates and mo- 
menta q a = S^Q/? and p a = S^P/?, we obtain a trun- 
cated form of the cavity Hamiltonian 

n cav = l^2(p 2 a +LJ 2 a q 2 a ), (30) 



where the new coordinate and momentum operators sat- 
isfy the appropriate commutation relations. Finally, the 
DHC modes are given by 



D(r,t) = V2Y1 

a, (3 

= V2^2u a F a (r)S ai 3 qp(0)cos(upt). 

a,(3 



(31) 



This means that the spatial distribution of the funda- 
mental mode j3 = 1 is 



D a (r) = V2j2^F a (r)S al . 



(32) 
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This expression represents of the DHC mode in terms of 
the bound PCW mode basis with the nRW wave terms 
retained. Using this representation the total electromag- 
netic energy in the fundamental cavity mode can be ap- 
proximated in terms of L a p and its eigenvectors as 



U> 



2 o2 

J rv°rv - 



(31 



L, 



a{3 



3(3a 



(33) 



For the cavities we have examined, we have found that 



the contribution of the second term in (33) is negligible. 
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FIG. 3: (Color online) Schematics of the two different cavity 
types we consider, (a) Photosensitive cavity: yellow shading 
represents a local change in the refractive index of the PCW 
slab creating the cavity, (b) Fluid infiltrated cavity: the red 
shading indicates a change in refractive index of the holes. 
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FIG. 2: (Color online) (a) TE-like band diagram of the PCW 
underlying Cavity 2. Broken blue line is the even PCW band 
used in the mode expansion. Red lines and shading show 
other modes, (b)-(e) |y • D&(r)| for different bound PCW 
modes indicated in (a) the crosses. The color scale is linear. 



B. Mode calculations 

In this section, we compare profiles of DHC modes 
computed using the Hamiltonian method with those com- 
puted using FDTD. Double- heterostructure cavities have 
been realized in different ways H 03 E3 H3 EMU ES]. 
Here we investigate the two geometries shown schemati- 
cally in Fig. [3j In the first (Fig. 3ja)), the photosensitive 
cavity, the refractive index of the PC slab is uniformly 
increased, often through a photo-induced refractive index 
change [TTJ [18] . The second is the fluid infiltrated cavity 
where the refractive index of the holes is increased, as 
can be achieved by fluid infiltration [20j [2TJ [35] . 

For these two cavity geometries we use three different 
sets of parameters for our calculations: 



1. A fluid infiltrated cavity based on a WO. 98 PCW 
with background index n& = 3.46 (consistent with 
silicon at A ~ 1500 nm, where A is the wavelength), 
air hole radius a = 0.26d, where d is the period, and 
slab thickness t = 0A9d. The cavity is introduced 
by increasing the refractive index of the holes by 
Ani = 0.2,0.4,0.6. 

2. A photosensitive cavity based on a Wl PCW with 
background index = 2.7 (the refractive index of 
some photosensitive chalcogenide glasses [36 J, air 
hole radius a = 0.3d and slab thickness t = 0.7d. 
The cavity is written by uniformly increasing the 
refractive index of the slab by An p = 0.02, 0.04. 

3. A photosensitive cavity based on a W0. 9 PCW with 
background index = 3.46, air hole radius a = 
0.3d and slab thickness t = 0.7 d. The cavity is 
introduced by uniformly increasing the refractive 
index of the slab by An p = 0.02, 0.04. 

Cavity 1 is similar to the experimental geometry in [2T] 
and Cavity 2 is similar to that in [18 . Geometries such 
as Cavity 3 have not been realized but could result from 
ion bombardment of a semiconductor [37j [38] . 

Figure [4] compares DHC modes computed using 
our Hamiltonian approach with those computed using 
FDTD. Figures 4ja),(d) and (g) show a z = slice of 
the electric field components E y computed by solving 
Eq. ( |27| ). Figures [4jb),(e),(h) are similar, but are along a 
y = slice. There is good agreement both in the modes' 
envelope and in the underlying rapid oscillations. Con- 
sidering the modes' complete 2D cross-section leads to 
the same conclusion. Figures gc),(f),(i) show the Bloch 
mode coefficient |5 C) fci| for the modes lying below the 
light line. In all cases the magnitude of the Bloch modes 
|£c 5 fci| is sufficiently small at the left edges of these plots 
to justify our truncation or Bloch mode basis to those 
below the light cone. However, Fig. [4^f) indicates that 
this approximation may break down if the cavity is made 
shorter than L = lOd. This is because the even PCW 
modes of slabs with a background index of = 2.7 have 
a higher frequency than those with rib = 3.46 so there 
are fewer non-radiative basis functions available for the 
mode expansion. 
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FIG. 4: Mode profiles of DHC modes taken at z = and their Bloch mode coefficients. Left column: E y field computed 
using the Hamiltonian formulation. Centre column: comparison of \E y \ (taken at y = z = 0) computed using the Hamiltonian 
formulation (broken red curve) and using FDTD (blue curve). Right column: magnitude of the Bloch mode coefficients (blue 
dots) l^cfcil for half of the BZ. The continuous curve is an interpolation of |5c,/ei| and is included to aid the eye. The red 
shading shows the value of k where the PCW band is inside the light line, (a)-(c) Cavity 1 with length L — 12d and index 
change Am = 0.2; (d)-(f) Cavity 2 with length L = lOd and An p = 0.02; (g)-(i) Cavity 3 with length L = 8d and An p = 0.02. 



Closer inspection of Fig. [4^b) shows a slight discrep- 
ancy between the modal widths from our theory and from 
FDTD. This is surprising since we expect the theory to 
be well-suited for modes such with a full-width at half 
maximum of approximately lOd and a highly localized 
Bloch mode composition (Figure |4jc)). However, in con- 
trast to photosensitive cavities, fluid infiltrated cavities 
are formed by large refractive index changes (An^ =0.2 
vs An p = 0.02) in regions where the electric fields of the 
PCW modes are weak. The weakness stems from the di- 
electric nature of the PCW modes, so they are localised 
in the slab rather than the holes, and from the absence of 
holes in the waveguide region where the field is strongest. 
We believe that our results for the fluid infiltrated cav- 
ity could be improved by including more basis modes in 
the mode expansion, for example modes which are not 
bound to the PCW, or higher order PCW modes. How- 
ever, our results are sufficiently accurate for qualitative 
insight into the modes and their radiation patterns. 



III. RADIATION FIELD OF DHC MODES 



Although D a (r) provides a good approximation to the 
field of the cavity mode, it is non-radiating. This implies 
that the Fourier transform of D a (r) has no Fourier com- 
ponents within the light cone. On the other hand, we 
can define a polarization field using D a (r), 



P a (r) =T(r)D a (r), 



(34) 



where 



r(r) 

f(r) 



<r) - 1 
e(r) ' 

e(r) - 1 
e(r) ' 



(35) 



We now have a good approximation for the shape of the 
cavity mode D a (r) and its frequency cDi, but the mode 
does not radiate. In the next section, these quantities are 
used as ingredients to formulate a perturbative treatment 
for computing the polarization field of DHC modes inside 
the light cone. These polarization fields are then used to 
compute the Q factor and far- field radiation patterns. 



with T(r) = r(r) + f(r), which defines f(r). The po- 
larization field P a (r) consists of a part T(r)D a (r) which 
does not radiate because D a (r) is entirely outside the 
light cone and T(r) has the periodicity of the PCW, and a 
radiating part T(r)D a (r). However, T(r) is not periodic 
and thus when multiplied by D a (r) does have Fourier 
components inside the light cone. 
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A. Green tensor for calculating radiation 

Before continuing we briefly review a formalism intro- 
duced earlier [29] [39] [40] for computing radiation fields. 
All material properties are placed in a polarization term, 
and the dynamic Maxwell equations with harmonic time 
dependence read 



iuB(r) = V x E(r), 
-^E(r) = V x B(r) + ^/i P(r). 



(36) 



Here the polarization field is a source term and we wish to 
find E(r) and B(r) for a given polarization distribution 
P(r). This implies that we require a Green tensor that 
propagates a field emitted by a polarization source oscil- 
lating at frequency uj. We adopt an earlier formulation 
[29] that expresses the Green tensor using the variables 
(k x ,k y ,z). In terms of the cavity mode, this distinguishes 
between modes that radiate k^ + ky < /cq, and those that 
are bound to the slab k\ + ky > /cq, where ko = uj/c. 

Taking k, = (k xi k y ) and R = (x,y), the relationship 
between the E(r) and P(r) is 



E(r) 



dz'Gfaz- z')>P(k,z'), (37) 



where P(/s, z) is the Fourier transform of P(r) in the x 
and y variables, and [29] 

G(k, z -z';lo) = ^ (e iw ^-^e(z - *')(§§ + p+p+) + 



-iw(z — z') 



Here, 0(z) is the Heaviside-step function, w = \/^o~ — n * 
where n = and s and p± are the unit vectors for the 
conventional s and p polarised plane waves respectively, 
where + and — refer to upward and downward propaga- 
tion. The definition of w includes the condition that \[X 
is defined so that ImV^ > 0, and if ImV^ = 0, then we 
take Yle\fX > 0. The directions of s and p± depend on 
k and are related to cartesian components by 



(39) 



where 0± points in the direction of propagation. These 
vectors form the orthogonal triads (p + ,s,z> + ) and 
(p_, s, z>_) for upward and downward propagating plane 
waves respectively. This is shown schematically for up- 
ward propagating modes in Fig. [5] The Green tensor 
contains the outgoing wave condition as required. The 
fields for a structure of finite thickness are obtained by 
convolving the Green tensor with the polarization distri- 
bution in the z-direction for each value of n. The spatial 
distribution of fields above and below the slab is then 
obtained by inverse Fourier transform. 
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FIG. 5: (Color online) Relationship between the cartesian 
wavevector components k = k x it + k y y + wz and its s and 
p polarised components for upward propagating plane waves. 
For + ky > ko the diagram is a schematic, since i>+, p+, 
and k are all complex. 



Returning to computing the radiation properties, we 
take the PC slab with thickness t to lie in the region 
between —t/2 < z < t/2. Using the Green tensor (38), 
the electric field above the slab is 



E+(r) 



idn 

2-KW 



e+(«), 



where 



e+(n) = se s + (K) + p + e^(K), 



6{z' - z){ss + p_p-)) - j- o 5(z - z')zz. (38) and 



e' + (n) 



k ° s- / dz'e- iwz 'P{K,z') 



e p + (n) = ^- P+ • J dz'e- iwz 'P(K,z'), 



(40) 



(41) 



(42) 



with an equivalent expression below the slab. The energy 
carried away from the cavity at any height above the 
slab is fixed. Furthermore, an expression for the field at 
any plane above the slab, say z = z$, can be used to 
propagate it to any other value of z > z$. This can be 
used to reproduce the results of Englund et al [25] . 

The asymptotic expression for the far-field in spherical 
polar coordinates is 



E(r,0,0)~e+(*)- 



ik r 



(43) 



where R, = k^r - (xx + yy) = &o (sin cos 0, sin sin (j)) , 
where (0, (j)) are the declination and azimuthal angles re- 
spectively. Here, r is a unit vector that identifies the 
direction in which we let r — >• oo, and therefore k is the 
projection of r onto the xy plane multiplied by fco- This 
then shows that each value of (k X: k y ) inside the light 
cone corresponds to radiation in a particular direction 
(0,(j)). The time averaged far- field Poynting vector for 
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an upward travelling wave in spherical coordinates is 

<S(r,M,t)> = [|e;(K)| 2 + |e^(K)| 2 ] u + , (44) 

where ( ) indicates a time average. The quality factor Q 
is the ratio of energy stored in the cavity mode and the 
energy lost per cycle and is given by 



The Green tensor is the same as that in Eq. (38), but 



U 



2 / d<t>dB sin 6S(6, 



(45) 



where S{6, (j>) = r 2 (S(r, 0, 0, t)) r is the radial component 
of the time-averaged Poynting vector as a function of 
the angles (0, (j>) and the integration is over the upper 
hemisphere. Here, U is the total energy of the field in 
the DHC mode, approximated by Eq. (33). Since the 



geometry is up-down symmetric with respect to the z- 
direction, the factor of 2 accounts for radiation emitted 
in the upper and lower hemisphere. 



B. The radiative polarization field of DHC modes 



Computing the far-field properties of DHC modes us- 



ing the theory presented in Section [Til A| requires an ex- 
pression for the Fourier components of the polarization 
field that lie inside the light cone. However, it would be 
wrong to approximate P(r) = P a (r) and to substitute 
its Fourier transform into Eqs. ( [42] ) . This is clear from 
P a (r) = f (r)D a (r) + f(r)D a (r). We consider the terms 
T(r), e(r) and D a (r) as being zeroth order, while the 
terms associated with the perturbation creating the cav- 
ity f(r) and e(r) are considered first-order small. Hence- 
forth, all parameters denoted by overbars are zeroth or- 
der small, while those with a tilde are first order small. 
Since f (r)D a (r) does not contribute to radiation, the 
radiative components of P a (r) are first-order small. Of 
course D a (r) is just an approximation for the actual field 
D(r) with corrections terms included by writing 



D(r) =D a (r)+D c (r), 



(46) 



where D c (r) contains first and higher order corrections. 
The polarization field is then given by 



P(r) =P°(r)+P c (r) 

= f (r)D°(r) + f(r)D a (r) + P c (r), 



(47) 



where the leading order term f (r)D a (r) is non-radiating 
and P c (r) contains first-order and higher order correc- 
tions. This means that P c (r) contains terms that are of 
the same order as T(r)D a (r), and thus an expression for 
P c (r) is required to compute the radiation fields. 

We begin by writing down a simple expression for a 
polarization field 



P(r) = e (e(r) - l)E(r) 



e (e(r)-l) J dv' G(v - v';w) ■ P(r'). 



(48) 



for brevity, we have expressed it using spatial variables. 
Since the cavity mode radiates and thus has a complex 
valued frequency, Eq. ( 48 ) has no solutions for real valued 



uj. However, from solving ([27]) we have a good approx- 
imation for the real part of the resonant frequency d>i, 
and similarly, P a (r) provides a starting point for approx- 
imating the polarization field. Using Eq. (47) we write 
Eq. (|48| as 



[P a (r) + P c (r)] = e (e(r) - 1) / dr f G(r - r', ^ + Q) 

.[P a (r')+P c (r')], (49) 

where uj is the complex first order correction to the fre- 
quency. We now take a Taylor expansion of the Green 
tensor in uj about Cj\ 

G(v-v'-u } ) = G(v-v'-C Jl ) + dG(v ~ V '' Q3l) C J + .... (50) 

UUJ 

To proceed, we write the Green tensor at each value of 
k as G(r — r'; uji) — G(r — r'; uj^) + G/c(r — r'; i.e. a 
sum of a Green tensor that propagates out fields with fre- 
quency cut and a correction term (which is a function of 
k). The frequency of DHC modes is typically close to the 
PCW band and thus (uji —0Jk)/&i <^ 1 for Bloch wavevec- 
tors that are below the light line. Therefore G/^r — r'; uji) 
is considered first-order small. The term 

e„(e(r) - 1) j dv'G{v - r';^) • f (r')D a (r') 
D°(r) 



= 6 (e(r)-l) 



e e(r) 



e (e(r) - 1) ^ c a f dr'G k (r - r'; oil) • f (r')F ce (r') ) 

(51) 



using Eq. (32) and c a = y/2(j a S a i. The first term on 
the RHS is 



T(r)e(r) 
e(r) 



(52) 



D°(r), 



and thus we have an equation for P c (r) 

P^) = ^D« (r) + ^«D«( r) 

+ eo(e(r) - 1) ^ c a J dr'G k (r - r'; oj,) ■ f (r')F Q (r') 



+ eo(e(r) -I) J dv'G{v - r'; wj) • [r(r')D a (r') + P c (r') 

/ / x ,\ - f , ,dG(r - r';o>i) 
+ e (e(r)-l)w / dv' 

f (r')D a (r') + f(r')D a (r') + P c (r')l + . . . . 

(53) 
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FIG. 6: (Color online) Schematic of the S filter in Fourier 
space. Red shading indicates the Fourier components that 
are kept, while white indicates Fourier components that are 
removed. The light cone is shown at the centre. 



We only take the terms that are first-order small in the 
tilde variables e(r), T(r), Gfc(r — r';u;), Co, and the first 
order polarization correction P c (r), giving an integral 
equation for the polarization field to first order 



PS(r) = ^fV(r) 



e (e(r) - 1) c a f dv'G k {r - r'; w x ) ■ f (r')F Q (r') 
e„(e(r) -I) j dv'G(y - r';^) • [f(r')D a (r') + PJ(r') 



+ e (e(r) - l)Q j dr' ^^^ ■ f (r')D»(r'). 



(54) 



To compute the radiative polarization field we only 
require terms with Fourier components inside the light 
cone, i.e. the first and third terms on the RHS of (54). 



It is tempting to simply apply a filter to both sides of 
Eq. (54) and eliminate all Fourier components outside 



the light cone. However, such a filter does not commute 
with multiplication by (e(r) — 1) and therefore can not be 
taken into the integral. This is because e(r) has the pe- 
riodicity of the PCW, and so mixes Fourier components 
separated by k x = 2iv/d. The requirements of a filter 
to commute with multiplication by e(r) are: it must be 
periodic in k x with period 2i\ / d, and it must be invariant 
under translations in k v . We therefore use 



§= f-^-e iK R y rect 
J (2tt) 2 ^ 



c\k x \ 
2cji 



• 27rm 



)/ 



dR'e 



-IK.R' 



(55) 

where rect(x) = 1 if \x\ < \ and is zero otherwise. This 
operates on a function by Fourier transforming, applying 
a filter in the Fourier domain and then inverse Fourier 
transforming. The Fourier filter is composed of an infinite 
series of rect functions with width 2Cj\ separated by k x = 
2ir/d as shown in Fig. [6j On application on both sides 
of Eq. ( |54| ), the filter removes all terms without Fourier 
components inside the light cone. This is because all 
terms without Fourier components inside the light cone 



also do not possess Fourier components separated from 
the light cone by k x = 2irm/d, where m is an integer. 
Equation (54) then becomes 



SP c 1 (t)=S j D a (r) + 

e (e(r) - 1) / dr'G(r - t , ;6j 1 )-S \f(/)B a (r') + P?(r 



(56) 



By defining Pj' rad (r) 



>c,rad 



SP?(r) ? as well as P r 1 ad (r) 



(r) + ST(r)D a (r), Eq. (56) becomes 



Pi ad (r) =S 
+ e (e(r)- 



r(r)e(r) 
e(r) 

) J dv' G(: 



+f(r) 



D a (r) 



(57) 



r-r^i) -P^^). 



This is a Fredholm equation of the second kind whose 
solution for P^ ad (r) gives a first order approximation for 
the radiative components of the DHC mode. The inho- 
mogeneous term evaluates to what degree the perturba- 
tion terms e(r) and T(r) couple Fourier components from 
the non-radiative approximation of the field D a (r) into 
the light cone. Once Eq. (57) is solved for P^ ad (r), the 
radiation can be obtained by computing the Poynting 
vector from Eqs. (40)- (45). Our approach for obtaining 
solutions to Eq. (57) is outlined in Appendix [Aj 



IV. RADIATION CALCULATION RESULTS 

In this section we present a comparison between Q 
factors and far-field radiation patterns obtained using the 
FAR and those computed using fully numerical FDTD 
calculations. 

We carried out our FDTD calculations using the com- 
mercial package Q- Finder (RSoft). The computational 
parameters were tailored to compute modes with differ- 
ent Q factors, but in all cases we used symmetry to re- 
duce the computation domain to 1/8 of the DHC. Our 
computation domain ranged from (x, z) = (0, 0, 0) to 

(x,y,z) = (25.5rf, 13 ^ 7 1.35(f), where d is the period. 
The spatial discretization ranged from (Ax, Ay, Az) = 
d/22 to (Ax, Ay, Az) = d/26, with adjustments to man- 
age the computation time. The temporal discretization 
was always set to cAt = Ax/2. Depending on the pa- 
rameters, these calculations took between 10-50 hours on 
a 32 core cluster for each data point in Fig. [7| 

In Fig. [7] we compare the Q factor for the three cavi- 
ties versus cavity length L computed using the FAR and 
FDTD. Once the basis functions for the underlying PCW 
have been computed, the Q factor calculations using the 
FAR typically takes less than 15 minutes per data point. 
For all three cavities we obtain good agreement in the 
trends of the Q factors versus cavity length between FAR 
and FDTD. For the theoretically calculated Q factors in 
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FIG. 7: (Color online) Quality factors versus cavity length L 
computed using our FAR method (red symbols) and those 
computed using FDTD (blue symbols), (a) Q factors for 
Cavity 1 with Am = 0.2,0.4,0.6. (b) Q factors for Cavity 
2 with Am = 0.02,0.04. (b) Q factors for Cavity 3 with 
Am = 0.02,0.04. 



Figure [7Jb) , we varied the length of the cavity in a con- 
tinuous manner, while we have only provided values com- 
puted using FDTD at a number of intervening points as 
it is impractical to compute all points using FDTD. Note 
the large oscillations in Q as the cavity length is varied. 
We return to this in Sect. El 

In Figure [7^b) , the difference between Q factors com- 
puted using FAR and FDTD is at most 30% (~ 2% for 
their logarithms), while in Figure [7^c) the discrepancy is 
at most 35% (~ 2.5% for their logarithms). On the other 
hand the agreement between FAR and FDTD for Cavity 
1 (Figure [T^a)) is not as impressive, and the discrepancy 
is at most a factor of 2 (6% for their logarithms). This is 
likely to be due to the fact that the field profiles computed 
using the bound mode basis have a slight discrepancy 
when compared with FDTD (see Figure Qb)). Since 
the Q factors here are very large, a small discrepancy 
in D a (r) may cause a large change in the radiation prop- 
erties. Nevertheless, the chief aim of our semi-analytic 
approach is to obtain qualitative information regarding 
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FIG. 8: (Color online) Far- field Poynting vectors (S r ) for Cav- 
ity 1 with Am =0.2 computed using the FAR (left column) 
and those computed using FDTD (right column), (a) Far- 
field radiation pattern for cavity length L = 4d, (b) L = 8d 
and (c) L = 12d. Colors as in Figure [2] Here (j) and 9 are the 
azimuthal and declination angles respectively. 



the general trends in Q factor as a function of cavity pa- 
rameters, and the results shown in Figure [7] indicate that 
the theory has achieved this goal. 

We now examine the far-field radiation patterns (the 
radial component of the Poynting vector S r ) for DHC 
modes. Figure [8] shows a comparison of far- field radia- 
tion patterns for Cavity 1 computed using the FAR (left 
column) and those using FDTD (right column). There is 
good agreement between the two sets of far-field patterns; 
both show that the number of lobes in the radiation pat- 
tern increases as the cavity becomes longer. The cavity 
with length L = Ad has particularly strong radiation in 
the vertical direction (6 ~ 0), which has been recently 
shown to be useful for exciting cavity modes from free 
In Section 



space 
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VI 



we provide designs for DHCs 
whose modes are engineered to emit vertically. 

Considering now the photosensitive cavity, the far-field 
radiation patterns of the modes of Cavity 2 are shown in 
Figure [9] Unlike the fluid infiltrated cavity, here the ra- 
diation pattern is predominantly directed towards large 
declination angles 0. The agreement between theory and 
FDTD is again good as both predict similar radiation di- 
rections and both provide the same trends for the number 
of lobes in the radiation pattern as the cavity length in- 
creases. The Q factors of Cavity 2 are larger than those 
in Cavity 1, even though the refractive index of Cavity 1 
is larger (715 = 3.46) than that of Cavity 2 (n& = 2.7) and 
the modes of Cavity 1 have envelope functions that vary 
more slowly than those of Cavity 2 (see Figure [4]). This is 



11 



(a) 1 
L = 6d 



(b) 

L = Sd 



(c) 1 
L = lOd 



run* 



Is" 
.s 

i -r 

_^i 

"en 

.3 

i -r 
i 

.3 



L.J 



(a) 1 
L = 6d 



~ l cos(0) sin(6>) 1 1 cos(0) sin(0) 1 



FIG. 9: (Color online) Far- field Poynting vectors (S r ) for Cav- 
ity 2 with An p = 0.02 computed using FAR (left column) and 
those computed using FDTD (right column), (a) Far- field ra- 
diation pattern for cavity length L = 6d, (b) L = 8d and (c) 
L = lOd. Colors as in Figure [2] 
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FIG. 10: (Color online) Far-field Poynting vectors (S r ) for 
Cavity 3 with An p = 0.02 computed using the FAR (left 
column) and those computed using FDTD (right column), 
(a) Far-field radiation pattern for cavity length L = 6d, (b) 
L — 8d and (c) L = lOd. Colors as in Figure [5] 



because changes in the background index couple the light 
much more weakly to the the light cone than changes in 
the refractive index of the holes. This is clear from ex- 
amining the Fourier components in the driving term in 
Eq. (57). We discuss this in more detail in Section fvj 

Finally, Figure [l0| shows the far-field radiation patterns 
computed for Cavity 3. Again there is good qualitative 
agreement between theoretical (left column) and numer- 
ical results (right column). The radiation patterns here 
resemble those computed for Cavity 2 in Figure [9j how- 
ever there are fewer lobes in the radiation pattern. This is 
because the frequency of the modes of Cavity 3 are lower 
than those for Cavity 2, and therefore the light cone oc- 
cupies a smaller region in Fourier space. Again, this can 
be observed through an examination of the Fourier com- 
ponents of the driving term in Eq. (57) within the light 



cone. We discuss this in the following Section. 



V. ANALYSIS OF THE DRIVING TERM 

Having established the quantitative capabilities of our 
theory, we now demonstrate the physical insight avail- 
able due to its analytic nature. In Eq. (57), the pa- 



rameters associated with the cavity geometry D a (r), 
e(r) and T(r), reside in the inhomogeneous driving term 



r(r)e(r) 
e(r) 



+r(r) 



(r) = A(r)D a (r), while the Green ten- 
sor ensures a self consistent interaction between dipoles. 
We now show that the far-field of the DHC modes 



can be understood from the Fourier components of the 
A(r)D a (r) term below the light line. 

We first examine the Q factor oscillations Fig. ff\b) 
associated with sub-period changes in cavity length. Ap- 
parently, the magnitude of the radiation is strongly af- 
fected by how the cavity perturbation cuts across the 
PCW modes. To analyze this we normalise the inhomo- 
geneous term A(r)D a (r) by dividing it by the square root 
of the total energy in the cavity mode y/U. The Fourier 
components of the term A(r)D a (r)/V^ inside the light 
cone then indicate the strength of the radiating polariza- 
tion field with respect to the total energy in the cavity 
mode. Figure [Tit a) snows a z = slice of the y compo- 
nent of A(r)D a (r)/V^ for a cavity with refractive index 
change An p = 0.02 and length L = Ad. Figure 11 'c) is 
similar, but for a longer cavity of length L = 4.8d. While 
the results seem similar, the Fourier components that are 
inside the light cone, shown in Figures [TT]fo) and (d), dif- 
fer strongly: the magnitude of the field inside the light 
cone is considerably larger for the longer cavity, and we 
therefore expect that this cavity radiates more strongly 
than the shorter cavity. This is confirmed in Fig. [7^b), 
which shows that the Q factor of the cavity with length 
L = Ad is approximately 8 times larger than that with 
length L = 4.8d. 

Figures 12'a)-(c) show z = slices of the y • D a (r) 
component of the inhomogeneous term for the three cav- 
ity types. The driving term for Cavity 1 shown in Fig- 
ure [l2[ a) has a Fourier transform (Figure |T2[d)) with a 
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FIG. 11: (Color online) Contour plot of the y component 
of the inhomogeneous term (z — slice) normalised by the 
total electromagnetic energy of the mode for Cavity 2 with 
An p = 0.02. (a) Cavity length L = Ad. (b) Cavity length 
L — A.Sd. (c) and (d) show the light cone components of the 
Fourier transforms of (a) and (b) respectively. 
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FIG. 12: (Color online) Contour plot of the y component of 
the inhomogeneous term [z — slice) for (a) Cavity 1 with 
Am = 0.2 and L = 4d, (b) Cavity 2 with An p = 0.02 and L = 
Ad, and (c) Cavity 3 with An p = 0.02 and L — Ad (d),(e),(f) 
show the light cone components of the Fourier transforms of 
(a),(b),(c) respectively. Color bars are linear scales. 



strong DC component, indicating strong radiation in the 
vertical direction consistent with the computed Poynt- 
ing vector (Figure [8^ a)). The strong vertical radiation 
means that this cavity design typically has a smaller Q 
factor than the photosensitive design. The driving term 
for Cavity 2 has a Fourier transform that is strongly 
peaked at the edges of the light cone (Fig. 12 'e)) and 



therefore the radiation is strongest at large declination 
angles, consistent with Fig. [9] There is also a subtle dif- 
ference between the far-fields for the two different pho- 
tosensitive cavities, i.e. for Cavities 2 and 3. The latter 
have a higher background index implying that its modes 
have lower frequencies, and consequently the light cone 
is smaller. This means that fewer features of the driving 
term overlap the light cone, explaining why Cavity 3 has 
far-field radiation patterns with fewer lobes (Figure 10) 
than Cavity 2 (Figure [9). 



VI. VERTICAL EMISSION 



Recent interest in engineering the radiation pattern of 
PC cavities [8j|4TJ[42] nas focused on cavities which emit 
radiation predominantly in the vertical. These not only 
enable the collection of light exiting the cavity, but al- 
low the cavity mode to be excited from free-space. These 
cavities were recently used in cavity QED [8 , 43 and har- 
monic generation experiments [14 . We show here that 
using the FAR we easily arrive at such designs for DHCs. 

As discussed in Sect. [V| the fundamental mode of a 
DHC to radiate predominantly in the vertical direction 
if the refractive index is such that A(r)D a (r) has a non- 
zero DC Fourier component. While we showed that this 
is so in a fluid-infiltrated cavity with L = Ad, this idea ap- 
plies more generally. Vertical radiation can be achieved 
by manipulating the holes in a different way. Figure [l3^a) 
shows a schematic of a fluid infiltrated cavity of length Ad 
with its radiation pattern computed using FDTD in Fig. 
[l3]fo). A schematic of a design where the hole radius is 
decreased is shown in Fig. [l3fc); the associated radia- 
tion pattern in Fig. 13 ^d), for a structure in which the 
hole radius was decreased from 0.27<i to 0.25d in a VF0.94 
silicon waveguide with thickness t = 0.45d, confirms the 
predominantly vertical emission. We also computed the 
radiation pattern for a 1^0.94 PCW with holes radius 
0.27d and thickness t = 0.45d, where the holes shown in 
Fig. [H|e) are shifted to that of a W0.98 PCW and the 
holes drawn with the dashed red lines are shifted a further 
0.02^^, i.e. to where the equivalent holes of a 1^1.02 
PCW would be. The radiation pattern in Fig. 13 'f ). All 



three designs have more than 70% of their radiated power 
within a declination angle of 30° (white circles in Figs. 
[l3^b),(d),(f)). The computed Q factors for the parame- 
ters in Figs 13'd),(f) ranged between 2 x 10 5 and 4 x 10 5 . 



Unlike previous designs based on L3 cavities [41], the Q 
factor of these designs can be controlled independently of 
the radiation pattern. The theoretical Q factor of these 
cavities can be increased by reducing the strength of the 
perturbation that creates the cavity, i.e. by decreasing 
the change in radius or hole shift. These cavity designs 
may have improve the performance of cavity-based ex- 
periments in harmonic generation and cavity QED. 
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FIG. 13: (Color online) Schematics of DHC designs that 
maximize vertical emission and their corresponding radiation 
patterns (S r )> (a)-(b) Fluid infiltrated cavity, (c)-(d) Cavity 
created by radius change, (e)-(f) Cavity created by hole shift. 
The holes drawn with dashed red lines are shifted more than 
others. The white circles in the radiation pattern corresponds 
to a declination angle of 30°. The holes shifts and radius 
changes in (c) and (e) are not to scale. Colors as in Figure [5] 



VII. DISCUSSION AND CONCLUSION 

Previously [24j [44] , DHC modes have been character- 
ized by an envelope function that modulates a rapidly 
varying Bloch mode. Here, we show that this picture is 
insufficient for explaining some of the physics underly- 
ing the results presented here, and that a multiple Bloch 
mode approach is required. 

We first compare the modal fields produced using an 
envelope function formulation with the fields obtained by 
solving Eq. (27). The key difference between the work 
here and the envelope function picture is that here we 
construct the cavity mode by superposing multiple Bloch 
modes, while in the envelope function approach the D 
field of the DHC mode has the form 



De„v(r) = f(x)O v/d (r), 



(58) 



where f(x) is the envelope function and D 7r / rf (r) is the 
displacement field of the band-edge Bloch mode. Figure 
14 'a) shows the electric field at a z = slice calculated 



for Cavity 2 with L = lOd and An p = 0.02 using our 



Hamiltonian formulation, while Figure 14 'b) shows this 
mode computed using the envelope function-based the- 
ory. Even though these modes have almost the same 
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FIG. 14: (Color online) Electric field E y (r) for Cavity 2 
with length L = lOd and An p = 0.02 computed using (a) 
the Hamiltonian formulation and (b) the envelope function- 
based theory, (c) A y = z — slice of the magnitude 
of the ^-component of the electric field \E y \ computed us- 
ing the Hamiltonian formulation (dashed red curve) and the 
envelope-based theory (green curve). The dashed lines show 
the physical length of the cavity 



characteristic length, as shown in Figure 14 'c), there is 



a clear difference. The mode computed using the Hamil- 
tonian has a chevron-like feature which is absent from 
the envelope-based calculation. This chevron-like feature 
arises because of the difference between the Bloch modes 



in the mode superposition in Eq. (32). The fields of the 



different PCW Bloch modes are illustrated in Figures 
|2jb)-(e) for the PCW underlying Cavity 2. The impor- 
tant feature here is that the Bloch modes at different 
Bloch wavevectors have different functional forms with 
respect to the y variable. This leads to the chevron-like 
feature of the field profile when these modes superposed 
to construct the cavity mode through Eq. (22). This is 



not accounted for in the single Bloch mode theory since 
the functional form of the DHC mode in Eq. (58) only 
contains a single Bloch mode D 7r / rf (r) modulated by an 
envelope function f{x) that is only a function of x. The 
^/-dependence of the cavity mode therefore only lies the 
band-edge Bloch mode D 7r /^(r). 

In the envelope function theory, changes in the length 
of the cavity only manifest themselves through the en- 
velope function f(x). Consequently, if D env (r) is used 
to compute the Fourier components of A(r)D env (r)/V^ 
in the light cone, changes in the cavity parameters only 
affect the k x distribution and therefore increases in cav- 
ity length would lead to a monotonic increase in the Q 
factor. This is because the characteristic length of f(x) 
increases as the cavity becomes longer, and hence the 



14 



Fourier transform of f(x) becomes narrower. The en- 
velope function theory cannot predict the differences in 
the Fourier components in Figs 11 'b) and (d) as these 
functions differ in both k x and k y . The Q oscillations 
are inherently caused by the interference of the different 
Bloch components of D a (r) in the Fourier components of 
the product A(r)D a (r)/V^ inside the light cone. 

The oscillations in Q, occurring as a result of superpo- 
sition of multiple Bloch modes, have a natural interpreta- 
tion in terms of Fabry- Perot resonances in the waveguide 
cavity: as the length of the cavity increases the wave- 
guide impedance at the end 'facets' (the plane where the 
perturbation ends) changes periodically. This leads to a 
change in the spacing of the fringes in the x-direction of 
Fourier space, which in turn leads to a periodic modula- 
tion of the Q factor. This supports the interpretation of 
Sauvan et al. [27 j, in which a dominant contribution to 
the loss of a DHC cavity arises from Fabry-Perot reflec- 
tions at the cavity boundaries. 

In conclusion, we have presented the first-principles 
FAR for calculating the near-field and far-field proper- 
ties of double-heterostructure cavity modes. Our theory 
is successful on two levels: it enables accurate numerical 
calculations of the Q factor and far-field radiation pattern 
of DHC modes, but more importantly, significant quali- 
tative insight into the far-field properties of DHC modes 
is gained by examining the inhomogeneous driving term 
in the integral equation (57). This theory has the capa- 



bility to not only speed up numerical calculations but, 
since it provides a direct link between a cavity geometry 
and its far-field properties, to greatly enhance the ability 
to design cavities with tailored far-field properties. We 
have shown this by providing designs for ultrahigh Q cav- 
ities whose radiation pattern has been engineered to emit 
predominantly in the vertical direction. 
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Appendix A: Solving the far-field polarization 
integral equation 



Here we present our strategy for solving Eq. (57) for 
the polarization field P^ ad (r) using a spatially discrete 
basis. Since D a (r) is constructed by superposing Bloch 
modes, the grid spacing is chosen such that it resolves 
most of the Fourier components of the Bloch modes from 
which it is composed. We found that a discretization 



of (Ax,Ay,Az) = (d/24, (V3d/2)/16, d/24), where d is 
the period, was sufficient for this purpose. The Bloch 
modes of the even PCW band were discretized with 
Akd = 27r0.02, which implies a computation domain 
length of bOd in the x-direction, i.e. the x-coordinate 
spans [— 25d, 25d]. In the ^/-direction our computation 
domain ranges between [— 16y/3d/2, 16y/3d/2], while in 
the z-direction we only require points within the slab, 
i.e. between [— 1/2, t/2], with t the slab thickness. Upon 
discretization, Eq. ( p7| ) becomes an inhomogeneous ma- 
trix equation 



(I - EG)p = a, 



(Al) 



where I is the identity operator, P^ ad (r) — > p, eo(e(r) — 
1) — >• E, J dr' G(r — r';a)i) — >• G and a contains all in- 
homogeneous terms. Therefore when thus represented, 
operators become matrices, vector fields become vectors, 
and a convolution with the Green tensor is represented 
as a matrix multiplication. In practice we compute the 
convolution in the x and y coordinates through multi- 
plications in Fourier space. Computing the convolution 
directly in the z-direction is feasible because the slab is 
thin. This implies that the matrix vector product Gp 
computes a disc rete version of Eqs. (37) and ( |38| ). 

To solve Eq. (Al) we need to compute (I — EG) _1 a. 
The difficulty of this inversion depends on the nature of 
the matrix (i.e. its symmetry properties and sparsity, 
etc.) and its size. Although the matrix is sparse, from 
our discretization and domain size given above, the vector 
a has approximately N a ~ 10 7 elements in each of its 
three components and therefore, to solve the problem, 
we would be required to invert a sparse matrix of size 
37V a x 37V a . 



Equation (57) has the form of the well-known discrete 

These systems are 



dipole scattering problem 
typically too large to be solved directly and instead iter- 
ative methods are used. When the size of the problem 
becomes too large either the iteration does not converge, 
or the number of iterations becomes so large such that 
computations become impractical, even with the most so- 
phisticated iterative method. No known iterative method 
can solve a system of equations with dimensions of 3 x 10 7 
within a reasonable computation time, so the size of the 
problem must be reduced. 

Chamet and Rahmani [49] recently tested different 
iterative methods for the discrete dipole problem in 
which a field scattering off a sphere with an electric 
and magnetic response. They showed that for a sphere 
that is discretized into N ~ 200000 points, leading to 
67V x 6 N) matrices, the iterative Generalized Product- 
type Bi- Conjugate Gradient method (GPBiCG) [5Ql [51] 
performed best in terms efficiency and robustness. Like 
most iterative solvers, this method works by minimiz- 
ing a residual, which defines how well a solution satisfies 
the matrix equation. For a trial solution of the matrix 
equation Ax = b, the residual is defined as = b — Ax^. 
The quality of the solution increases as ||r^|| — >• 0. Itera- 
tive methods also do not have the storage of matrix G, 
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as only matrix vector products need to be computed. 
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FIG. 15: (Color online) The logarithm of the residual versus 
iteration number for a GPBiCG algorithm applied to solve 
Eq. (|Al) for Cavity 2 with An p = 0.02 and length L = Ad. 



Since we are ultimately interested in solving for Fourier 
components of p inside the light cone we only require 
the slowly varying components of p, thereby reducing 
the size of our problem to that in [49]. Although the 
rapidly varying Fourier components of the field are re- 
quired when finding the DHC mode D a (r), the radiative 
Fourier components are inherently slowly varying, and 
therefore, for this part of the problem we do not need a 
fine discretization in the x — y directions. We reduced the 
discretization to (Ax, Ay, Az) = (d/4, (>/3d/2)/4, d/24), 
where Az is unchanged. With this discretization each 
matrix vector product takes approximately 5 seconds to 
compute on our M ATLAB code. We note that the bulk 
of the computation time in the GPBiCG is taken up by 
the two matrix vector products in each iteration. Fig- 
ure HU shows the value of the residual versus iteration 
number for Cavity 2 with a cavity length L = Ad and 
An = 0.02. This shows that even after 2000 iterations 
there is no sign of convergence. 

Since we cannot achieve convergence we choose to solve 
Eq. ( Al ) under an approximation: we neglect all coupling 



to those outside t he l ight cone. This means multiplying 
both sides of Eq. (LA1J) by a F ourier filter F, that removes 



ky > fcg- Equation (Al) then becomes 



F(I - EG)p = Fa. 



(A2) 



We expect solutions to Eq. ( A2 ) to be approximate solu- 



tions to Eq. ( Al ). This is because we have observed that 
multiplying a vector that only has Fourier components 
inside the light cone by EG, results in a vector that is 
still dominated by its Fourier components inside the light 
cone, i.e. the light cone Fourier components are weakly 
coupled to Fourier components outside the light cone. 
Typical examples of the residual and the Q-factor ver- 



sus iteration number when solving Eq. (A2) using GP- 
BiCG are shown in Figure 16 We consider the problem 
to be solved when the residual is < 10 -5 which, in this 
example, is achieved in 52 iterations. Figure [16] shows 



that the calculated Q-factor also converges. We fu rther 
test our solution by substituting it into Eq. (| A1|) and 
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FIG. 16: (Color online) (a) The logarithm of the residual 
and (b) the Q-factor versus ite ratio n number for a GPBiCG 
algorithm applied to solve Eq. (A2). 



between Fourier components from inside the light cone 



checking if it satisfies this equation for the Fourier com- 
ponents inside the light cone-if we wish to see if p^ solves 
Eq. ( Al ) , we compute 1 1 p t — (EGpt + a) 1 1 / 1 1 p t 1 1 and check 
that it is of the same order as the residual (i.e ~ 10 -5 ). 
We found this to be so for all solutions presented. 
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